MODULE smm
  USE GLOBVAR
  USE NRUTIL
  USE PROBABILITY
  USE MATRIX
  !USE MPI
  USE random
  USE STATISTICS
  IMPLICIT NONE
  
    CONTAINS
   REAL(8) FUNCTION mom(thetain) 
   !SUBROUTINE mom(theta_ini,beta_f,check,N_int,sc)
        IMPLICIT NONE
        REAL(8), INTENT(IN):: thetain(:) !initial guess
        REAL(8):: Ls(N_inds,N_testt,3),Ep(N_inds,N_testt,3),as,as2,as3,qqm(N_inds,N_testt),qq(3)
        REAL(8):: L0s(N_inds,3),eta(N_inds,N_testt,3),e_L0(n_inds,3), e_eta(n_inds,3),mom_simu(N_mom),dif(N_mom,1)!, latent(N_inds*N_testl,N_test),ps(N_testl,N_test)
        
        ! parameters
        REAL(8):: beta_mu(n_x_l0,3),beta_eta(n_x_eta,3),var_l(3),etav(3),Lbar(N_levell),var_le(N_levell),deltal(N_levell),&
        Cbar(N_levelc),var_ce(N_levelc),deltac(N_levelc),Fbar(N_levelf),var_fe(N_levelf),deltaf(N_levelf),beta_l(3,3),beta_e(3,3),beta_y(2,3)  
        !gbar(N_levelg),var_ge(N_levelg),deltag(N_levelg),&
           
        INTEGER:: Ds(N_inds,N_testt),indext(n_testt,2), artls(n_inds,11,20),artcs(n_inds,N_levelc,11),artfs(n_inds,N_levelf,6)!,artgs(n_inds,N_levelg,3) 
        INTEGER:: i,j,k,jj,kk,m,t,tt,n,l,mmm
        REAL(8),ALLOCATABLE:: AVSE(:),cors(:,:)
        REAL(8):: inters_l(n_inds,n_levell,20,3),interls_l(n_inds,n_levell,3),inters_c(n_inds,n_levelc,20,3),intercs_c(n_inds,n_levelc,3),inters_f(n_inds,n_levelf,6,3),interfs_f(n_inds,n_levelf,3)
        INTEGER:: duls_l(n_inds,n_levell),dus_l(n_levell),inds_l(N_inds,n_levell),ducs_c(n_inds,n_levelc),dus_c(n_levelc),inds_c(N_inds,n_levelc),dufs_f(n_inds,n_levelf),dus_f(n_levelf),inds_f(N_inds,n_levelf) 

        
        !!!!!!!!!!!!!!!!!!Simulate Moments
        CALL SET_SEED(1438914432,727629,19890414,111111)
        ! language
        beta_mu(1:n_x_l0,1)=thetain(1:n_x_l0)                                              ! Initial Latent skill
        DO i=2,4
          beta_mu(i,1)=1.40d0**(thetain(i))/(10.0d0*(1.0d0+1.40d0**(thetain(i))))
        END DO
        beta_mu(5,1)=1.0d0-1.40d0**(thetain(5))
        
        var_l(1)=4.0d0*(1.4d0**thetain(n_x_l0+1)/(1.0d0+1.4d0**thetain(n_x_l0+1)))     ! variance of the shoch of initial latent skill
        Lbar=1.0d0
        Lbar(1)=0.0d0
        Lbar(3:N_levell) =thetain(n_x_l0+2:n_x_l0+N_levell-1)                ! mimimum requriement for each difficulty level
        var_le=1.0d0
        DO i=3,n_levell
          var_le(i)=1.4d0**thetain(n_x_l0+N_levell-1+i-2)     ! variance of the shock at each difficulty level
        END DO
        beta_eta(1:n_x_eta,1)=thetain(n_x_l0+2*N_levell-3+1:n_x_l0+2*N_levell+n_x_eta-3)           ! x for eta
        etav(1)=1.40d0**thetain(n_x_l0+2*N_levell+n_x_eta-2)                           ! variance of shock for eta
        
        DELTAl=1.0D0
        DO i=3,N_levell
            deltal(i)=6*(1.10d0**thetain(n_x_l0+2*N_levell+n_x_eta-2+i-2))/(1.0d0+1.10d0**thetain(n_x_l0+2*N_levell+n_x_eta-2+i-2))  ! delta for each level
        END DO
        tt=n_x_l0+3*N_levell+n_x_eta-4
        beta_y(1:2,1)=thetain(tt+1:tt+2)
        tt=tt+2
        ! cognitive 
        beta_mu(1:n_x_l0,2)=thetain(tt+1:tt+n_x_l0)                                              ! Initial Latent skill
        DO i=2,4
          beta_mu(i,2)=1.40d0**(thetain(tt+i))/(10.0d0*(1.0d0+1.40d0**(thetain(tt+i))))
        END DO
        beta_mu(5,2)=1.0d0-1.40d0**(thetain(tt+5))
        
        var_l(2)=4.0d0*(1.4d0**thetain(tt+n_x_l0+1)/(1.0d0+1.4d0**thetain(tt+n_x_l0+1)))     ! variance of the shoch of initial latent skill
        Cbar=1.0d0
        Cbar(2:N_levelc) =thetain(tt+n_x_l0+2:tt+n_x_l0+N_levelc)                          ! mimimum requriement for each difficulty level
        var_ce=1.0d0
        var_ce(2:N_levelc-1)=1.4d0**thetain(tt+1+n_x_l0+N_levelc:tt+n_x_l0+2*N_levelc-2)     ! variance of the shock at each difficulty level
        beta_eta(1:n_x_eta,2)=thetain(tt+n_x_l0+2*N_levelc-1:tt+n_x_l0+2*N_levelc+n_x_eta-2)           ! x for eta
        etav(2)=1.40d0**thetain(tt+n_x_l0+2*N_levelc+n_x_eta-1)                           ! variance of shock for eta
        deltac=1.0d0
        deltac(2)=4*(1.10d0**thetain(tt+n_x_l0+2*N_levelc+n_x_eta))/(1.0d0+1.10d0**thetain(tt+n_x_l0+2*N_levelc+n_x_eta))  ! delta for each level
        DO i=5,8
            deltac(i)=4*(1.10d0**thetain(tt+n_x_l0+2*N_levelc+n_x_eta+i-4))/(1.0d0+1.10d0**thetain(tt+n_x_l0+2*N_levelc+n_x_eta+i-4))  ! delta for each level
        END DO
        DO i=10,n_levelc-1
            deltac(i)=4*(1.10d0**thetain(tt+n_x_l0+2*N_levelc+n_x_eta+i-5))/(1.0d0+1.10d0**thetain(tt+n_x_l0+2*N_levelc+n_x_eta+i-5))  ! delta for each level
        END DO
        tt=tt+n_x_l0+3*N_levelc+n_x_eta-6
        beta_y(1:2,2)=thetain(tt+1:tt+2)
        tt=tt+2
        ! fine motor
        beta_mu(1:n_x_l0,3)=thetain(tt+1:tt+n_x_l0)                                              ! Initial Latent skill
        DO i=2,4
          beta_mu(i,3)=1.40d0**(thetain(tt+i))/(10.0d0*(1.0d0+1.40d0**(tt+thetain(i))))
        END DO
        beta_mu(5,3)=1.0d0-1.40d0**(thetain(tt+5))
        
        var_l(3)=4.0d0*(1.4d0**thetain(tt+n_x_l0+1)/(1.0d0+1.4d0**thetain(tt+n_x_l0+1)))     ! variance of the shoch of initial latent skill
        fbar=1.0d0
        fbar(2:N_levelf) =thetain(tt+n_x_l0+2:tt+n_x_l0+N_levelf)                          ! mimimum requriement for each difficulty level
        var_fe=1.0d0
        var_fe(2:N_levelf)=1.4d0**thetain(tt+1+n_x_l0+N_levelf:tt+n_x_l0+2*N_levelf-1)         ! variance of the shock at each difficulty level
        beta_eta(1:n_x_eta,3)=thetain(tt+n_x_l0+2*N_levelf:tt+n_x_l0+2*N_levelf+n_x_eta-1)           ! x for eta
        etav(3)=1.40d0**thetain(tt+n_x_l0+2*N_levelf+n_x_eta)                                    ! variance of shock for eta
        deltaf=1.0d0
        DO i=3,N_levelf
            IF (i==3) THEN
                 deltaf(i)=4.0D0 +thetain(tt+n_x_l0+2*N_levelf+n_x_eta+i-2)
            ELSE IF (i==6) THEN
                 deltaf(i)=4.0D0 +thetain(tt+n_x_l0+2*N_levelf+n_x_eta+i-2)
            ELSE
                 deltaf(i)=4*(1.40d0**thetain(tt+n_x_l0+2*N_levelf+n_x_eta+i-2))/(1.0d0+1.40d0**thetain(tt+n_x_l0+2*N_levelf+n_x_eta+i-2))  ! delta for each level
            END IF
        END DO
        
        tt=tt+n_x_l0+3*N_levelf+n_x_eta-2
        beta_y(1:2,3)=thetain(tt+1:tt+2)
        tt=tt+2
        DO i=1,3
            DO j=1,3
                beta_l(i,j)=thetain(tt+(j-1)*3+i)!2.0d0*(2.0d0*(1.1d0**(thetain(tt+(j-1)*3+i))/(1.0d0+1.1d0**(thetain(tt+(j-1)*3+i))))-1.0d0)
            END DO
        END DO
        tt=tt+9
        DO i=1,3
            DO j=1,3
                beta_e(i,j)=thetain(tt+(j-1)*3+i)!2.0d0*(2.0d0*(1.1d0**(thetain(tt+(j-1)*3+i))/(1.0d0+1.1d0**(thetain(tt+(j-1)*3+i))))-1.0d0)
            END DO
        END DO
        
        tt=tt+9
        !!!!!!!! ADD MEASUREMENT SHOCK!!!!!
        DO i=1,3
            qq(i)=1.40d0**(thetain(tt+i))/(1.0d0+1.40d0**(thetain(tt+i)))
        END DO
        

        
        
        ! generate shock terms
        
        DO i=1,N_inds
            DO j=1,3
              e_L0(i,j) =Sample_NORMAL(0.0d0,var_l(j))
              e_eta(i,j)=Sample_Normal(0.0d0,etav(j))
            END DO
        END DO
        
        ! GENERATE INITIAL LATENT VALUES AND 
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            DO j=1,3
              L0s(i,j)=LOG(1.10d0**(DOT_PRODUCT(L0(k,:),beta_mu(:,j))+e_L0(i,j))) ! NOTICE HERE IT IS LOG VALUES
            END DO
        END DO
        
        ! GENERATE MEASUREMENT ERRORS
        DO i=1,N_inds
            Do j=1,N_testt
                qqm(i,j)=Sample_Uniform(0.0d0,1.0d0)
            END DO
        END DO

        
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
       ! SIMULATE LATENT SKILL PROCESS ! HERE ONLY FOR LANGUAGE SKILLS 
        
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            DO j=1,N_testt
                IF (j<firstt(k,1)) THEN
                    Ls(i,j,1:3)=-99.0d0
                ELSE IF (j==firstt(k,1) .and. m_t(k,j)==1) THEN
                    DO m=1,3
                        IF (x_testt(k,j,1)<-90) THEN
                            eta(i,j,m)=1.40d0**(e_eta(i,m))/((1.0d0+1.40d0**(e_eta(i,m))))
                        ELSE
                            eta(i,j,m)=1.40d0**(DOT_PRODUCT(x_testt(k,j,:),beta_eta(:,m))+e_eta(i,m))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testt(k,j,:),beta_eta(:,m))+e_eta(i,m))))
                        END IF
                    END DO
                    IF (levell(j)>-90 .and. levelc(j)<-90) THEN
                      kk=levell(j)
123                   as=Sample_Normal(0.0d0,var_le(kk))
                      as2=Sample_Normal(0.0d0,var_le(kk))
                      as3=Sample_Normal(0.0d0,var_le(kk))
                      Ls(i,j,1)=(DOT_PRODUCT(L0s(i,:),beta_l(:,1)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+beta_y(1,1)*month(j)+beta_y(2,1)*DBLE(time(j))+as
                      Ls(i,j,2)=(DOT_PRODUCT(L0s(i,:),beta_l(:,2)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+beta_y(1,2)*month(j)+beta_y(2,2)*DBLE(time(j))+as2
                      Ls(i,j,3)=(DOT_PRODUCT(L0s(i,:),beta_l(:,3)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3

                      !IF (minval(Ls(i,j,1:4))<0.0d0) GO TO 123
                    ELSE IF (levell(j)<-90 .and. levelc(j)>-90) THEN
                        kk=levelc(j)
131                     as=Sample_Normal(0.0d0,var_ce(kk))
                        as2=Sample_Normal(0.0d0,var_ce(kk))
                        as3=Sample_Normal(0.0d0,var_ce(kk))
                        Ls(i,j,1)=(DOT_PRODUCT(L0s(i,:),beta_l(:,1)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+beta_y(1,1)*month(j)+beta_y(2,1)*DBLE(time(j))+as
                        Ls(i,j,2)=(DOT_PRODUCT(L0s(i,:),beta_l(:,2)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+beta_y(1,2)*month(j)+beta_y(2,2)*DBLE(time(j))+as2
                        Ls(i,j,3)=(DOT_PRODUCT(L0s(i,:),beta_l(:,3)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3
                        !IF (minval(Ls(i,j,1:4))<0.0d0) GO TO 131
                    ELSE IF (levell(j)>-90 .and. levelc(j)>-90) THEN
                        kk=levelc(j)
                        tt=levell(j)
                        as=Sample_Normal(0.0d0,var_ce(kk))
                        as2=Sample_Normal(0.0d0,var_le(tt))
                        as3=Sample_Normal(0.0d0,var_ce(kk))
                        Ls(i,j,1)=(DOT_PRODUCT(L0s(i,:),beta_l(:,1)))+deltal(tt)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+as2
                        Ls(i,j,2)=(DOT_PRODUCT(L0s(i,:),beta_l(:,2)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+as
                        Ls(i,j,3)=min((DOT_PRODUCT(L0s(i,:),beta_l(:,3)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+as3,(DOT_PRODUCT(L0s(i,:),beta_l(:,3)))+deltal(tt)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+as3)
                    ELSE IF (levell(j)<-90 .and. levelc(j)<-90 .and. levelf(j)>-90) THEN
                        kk=levelf(j)
139                     as =Sample_Normal(0.0d0,var_fe(kk))
                        as2=Sample_Normal(0.0d0,var_fe(kk))
                        as3=Sample_Normal(0.0d0,var_fe(kk))
                        Ls(i,j,1)=(DOT_PRODUCT(L0s(i,:),beta_l(:,1)))+deltaf(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+beta_y(1,1)*month(j)+beta_y(2,1)*DBLE(time(j))+as
                        Ls(i,j,2)=(DOT_PRODUCT(L0s(i,:),beta_l(:,2)))+deltaf(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+beta_y(1,2)*month(j)+beta_y(2,2)*DBLE(time(j))+as2
                        Ls(i,j,3)=(DOT_PRODUCT(L0s(i,:),beta_l(:,3)))+deltaf(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3
                        !IF (minval(Ls(i,j,1:4))<0.0d0) GO TO 139
                    ELSE IF (levell(j)<-90 .and. levelc(j)<-90 .and. levelf(j)<-90 .AND. levelg(j)>-90) THEN
                        Ls(i,j,1:3)=-99.0d0
                    END IF
                ELSE IF (j==firstt(k,1) .and. m_t(k,j)==-1) THEN
                    Ls(i,j,1:3)=-99.0d0
                ELSE IF (j>firstt(k,1) .and. Ls(i,j-1,1)<0.0d0 .and. m_t(k,j)==-1) THEN
                    Ls(i,j,1:3)=-99.0d0
                ELSE IF (j>firstt(k,1) .and. Ls(i,j-1,1)<0.0d0 .and. m_t(k,j)==1) THEN
                    DO m=1,3
                        IF (x_testt(k,j,1)<-90) THEN
                            eta(i,j,m)=1.40d0**(e_eta(i,m))/((1.0d0+1.40d0**(e_eta(i,m))))
                        ELSE
                            eta(i,j,m)=1.40d0**(DOT_PRODUCT(x_testt(k,j,:),beta_eta(:,m))+e_eta(i,m))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testt(k,j,:),beta_eta(:,m))+e_eta(i,m))))
                        END IF
                    END DO
                    IF (levell(j)>-90 .and. levelc(j)<-90) THEN
                      kk=levell(j)
164                   as=Sample_Normal(0.0d0,var_le(kk))
                      as2=Sample_Normal(0.0d0,var_le(kk))
                      as3=Sample_Normal(0.0d0,var_le(kk))
                      Ls(i,j,1)=(DOT_PRODUCT(L0s(i,:),beta_l(:,1)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+beta_y(1,1)*month(j)+beta_y(2,1)*DBLE(time(j))+as
                      Ls(i,j,2)=(DOT_PRODUCT(L0s(i,:),beta_l(:,2)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+beta_y(1,2)*month(j)+beta_y(2,2)*DBLE(time(j))+as2
                      Ls(i,j,3)=(DOT_PRODUCT(L0s(i,:),beta_l(:,3)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3
                      !IF (minval(Ls(i,j,1:4))<0.0d0) GO TO 164
                    ELSE IF (levell(j)<-90 .and. levelc(j)>-90) THEN
                        kk=levelc(j)
172                     as=Sample_Normal(0.0d0,var_ce(kk))
                        as2=Sample_Normal(0.0d0,var_ce(kk))
                        as3=Sample_Normal(0.0d0,var_ce(kk))
                        Ls(i,j,1)=(DOT_PRODUCT(L0s(i,:),beta_l(:,1)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+beta_y(1,1)*month(j)+beta_y(2,1)*DBLE(time(j))+as
                        Ls(i,j,2)=(DOT_PRODUCT(L0s(i,:),beta_l(:,2)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+beta_y(1,2)*month(j)+beta_y(2,2)*DBLE(time(j))+as2
                        Ls(i,j,3)=(DOT_PRODUCT(L0s(i,:),beta_l(:,3)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3
                        !IF (minval(Ls(i,j,1:4))<0.0d0) GO TO 172
                    ELSE IF (levell(j)>-90 .and. levelc(j)>-90) THEN
                        kk=levelc(j)
                        tt=levell(j)
                        as=Sample_Normal(0.0d0,var_ce(kk))
                        as2=Sample_Normal(0.0d0,var_le(tt))
                        as3=Sample_Normal(0.0d0,var_ce(kk))
                        Ls(i,j,1)=(DOT_PRODUCT(L0s(i,:),beta_l(:,1)))+deltal(tt)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+beta_y(1,1)*month(j)+beta_y(2,1)*DBLE(time(j))+as2
                        Ls(i,j,2)=(DOT_PRODUCT(L0s(i,:),beta_l(:,2)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+beta_y(1,2)*month(j)+beta_y(2,2)*DBLE(time(j))+as
                        Ls(i,j,3)=min((DOT_PRODUCT(L0s(i,:),beta_l(:,3)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3,&
                            (DOT_PRODUCT(L0s(i,:),beta_l(:,3)))+deltal(tt)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3)
                        
                    ELSE IF (levell(j)<-90 .and. levelc(j)<-90 .and. levelf(j)>-90) THEN
                        kk=levelf(j)
180                     as=Sample_Normal(0.0d0,var_fe(kk))
                        as2=Sample_Normal(0.0d0,var_fe(kk))
                        as3=Sample_Normal(0.0d0,var_fe(kk))
                        Ls(i,j,1)=(DOT_PRODUCT(L0s(i,:),beta_l(:,1)))+deltaf(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+beta_y(1,1)*month(j)+beta_y(2,1)*DBLE(time(j))+as
                        Ls(i,j,2)=(DOT_PRODUCT(L0s(i,:),beta_l(:,2)))+deltaf(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+beta_y(1,2)*month(j)+beta_y(2,2)*DBLE(time(j))+as2
                        Ls(i,j,3)=(DOT_PRODUCT(L0s(i,:),beta_l(:,3)))+deltaf(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3

                        !IF (minval(Ls(i,j,1:4))<0.0d0) GO TO 180
                    ELSE IF (levell(j)<-90 .and. levelc(j)<-90 .and. levelf(j)<-90 .AND. levelg(j)>-90) THEN
                        Ls(i,j,1:3)=-99.0d0
                    END IF
                    
                ELSE IF (j>firstt(k,1) .and. Ls(i,j-1,1)>=0.0d0 .and. m_t(k,j)==-1) THEN
                    IF (levell(j)>-90 .and. levelc(j)<-90) THEN
                        Ls(i,j,1)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,1)))
                        Ls(i,j,2)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,2)))
                        Ls(i,j,3)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,3)))
                    ELSE IF (levell(j)<-90 .and. levelc(j)>-90) THEN
                        Ls(i,j,1)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,1)))
                        Ls(i,j,2)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,2)))
                        Ls(i,j,3)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,3)))
                    ELSE IF (levell(j)>-90 .and. levelc(j)>-90) THEN
                        Ls(i,j,1)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,1)))
                        Ls(i,j,2)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,2)))
                        Ls(i,j,3)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,3)))
                    ELSE IF (levell(j)<-90 .and. levelc(j)<-90 .and. levelf(j)>-90) THEN
                        Ls(i,j,1)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,1)))
                        Ls(i,j,2)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,2)))
                        Ls(i,j,3)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,3)))
                    ELSE IF (levell(j)<-90 .and. levelc(j)<-90 .and. levelf(j)<-90 .AND. levelg(j)>-90) THEN
                        Ls(i,j,1)=Ls(i,j-1,1)
                        Ls(i,j,2)=Ls(i,j-1,2)
                        Ls(i,j,3)=Ls(i,j-1,3)
                    END IF
                ELSE IF (j>firstt(k,1) .and. Ls(i,j-1,1)>=0.0d0 .and. m_t(k,j)==1) THEN
                    DO m=1,3
                        IF (x_testt(k,j,1)<-90) THEN
                            eta(i,j,m)=1.40d0**(e_eta(i,m))/((1.0d0+1.40d0**(e_eta(i,m))))
                        ELSE
                            eta(i,j,m)=1.40d0**(DOT_PRODUCT(x_testt(k,j,:),beta_eta(:,m))+e_eta(i,m))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testt(k,j,:),beta_eta(:,m))+e_eta(i,m))))
                        END IF
                    END DO
                    IF (levell(j)>-90 .and. levelc(j)<-90) THEN
                      kk=levell(j)
204                   as=Sample_Normal(0.0d0,var_le(kk))
                      as2=Sample_Normal(0.0d0,var_le(kk))
                      as3=Sample_Normal(0.0d0,var_le(kk))
                      Ls(i,j,1)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,1)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+beta_y(1,1)*month(j)+beta_y(2,1)*DBLE(time(j))+as
                      Ls(i,j,2)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,2)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+beta_y(1,2)*month(j)+beta_y(2,2)*DBLE(time(j))+as2
                      Ls(i,j,3)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,3)))+deltal(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3

                      !IF (minval(Ls(i,j,1:4))<0.0d0) GO TO 204
                    ELSE IF (levell(j)<-90 .and. levelc(j)>-90) THEN
                        kk=levelc(j)
212                     as=Sample_Normal(0.0d0,var_ce(kk))
                        as2=Sample_Normal(0.0d0,var_ce(kk))
                        as3=Sample_Normal(0.0d0,var_ce(kk))
                        Ls(i,j,1)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,1)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+beta_y(1,1)*month(j)+beta_y(2,1)*DBLE(time(j))+as
                        Ls(i,j,2)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,2)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+beta_y(1,2)*month(j)+beta_y(2,2)*DBLE(time(j))+as2
                        Ls(i,j,3)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,3)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3


                        !IF (minval(Ls(i,j,1:4))<0.0d0) GO TO 212
                    ELSE IF (levell(j)>-90 .and. levelc(j)>-90) THEN
                        kk=levelc(j)
                        tt=levell(j)
                        as=Sample_Normal(0.0d0,var_ce(kk))
                        as2=Sample_Normal(0.0d0,var_le(tt))
                        as3=Sample_Normal(0.0d0,var_ce(kk))
                        Ls(i,j,1)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,1)))+deltal(tt)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+beta_y(1,1)*month(j)+beta_y(2,1)*DBLE(time(j))+as2
                        Ls(i,j,2)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,2)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+beta_y(1,2)*month(j)+beta_y(2,2)*DBLE(time(j))+as
                        Ls(i,j,3)=min((DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,3)))+deltac(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3,&
                            (DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,3)))+deltal(tt)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3)

                    ELSE IF (levell(j)<-90 .and. levelc(j)<-90 .and. levelf(j)>-90) THEN
                        kk=levelf(j)
220                     as=Sample_Normal(0.0d0,var_fe(kk))
                        as2=Sample_Normal(0.0d0,var_fe(kk))
                        as3=Sample_Normal(0.0d0,var_fe(kk))
                        Ls(i,j,1)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,1)))+deltaf(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,1))+beta_y(1,1)*month(j)+beta_y(2,1)*DBLE(time(j))+as
                        Ls(i,j,2)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,2)))+deltaf(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,2))+beta_y(1,2)*month(j)+beta_y(2,2)*DBLE(time(j))+as2
                        Ls(i,j,3)=(DOT_PRODUCT(Ls(i,j-1,:),beta_l(:,3)))+deltaf(kk)*DOT_PRODUCT(eta(i,j,:),beta_e(:,3))+beta_y(1,3)*month(j)+beta_y(2,3)*DBLE(time(j))+as3

                        !IF (minval(Ls(i,j,1:4))<0.0d0) GO TO 220
                    ELSE IF (levell(j)<-90 .and. levelc(j)<-90 .and. levelf(j)<-90 .AND. levelg(j)>-90) THEN
                        !kk=levelg(j)
!228                     !as=Sample_Normal(0.0d0,var_ge(kk-1))
                        Ls(i,j,1)=Ls(i,j-1,1)
                        Ls(i,j,2)=Ls(i,j-1,2)
                        Ls(i,j,3)=Ls(i,j-1,3)

                        !IF (minval(Ls(i,j,1:4))<0.0d0) GO TO 228
                    END IF
                END IF
            END DO
        END DO
        
        ! GENERATE TASK PERFORMANCE !
        ds=-99
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            DO j=1,N_testt
                IF (levell(j)>-90 .and. levelc(j)<-90) THEN
                   kk=levell(j)
                   IF (m_t(k,j)==-1) THEN
                       DS(i,j)=-99
                   ELSE IF (EXP(Ls(i,j,1))>=SUM(Lbar(1:kk)) .and. m_t(k,j)==1) THEN
                       IF (qqm(i,j)>=qq(1)) THEN
                           DS(i,j)=1
                       ELSE IF (qqm(i,j)<qq(1)) THEN
                           DS(i,j)=0
                       END IF
                       
                   ELSE IF (EXP(Ls(i,j,1))<SUM(Lbar(1:kk)) .and. m_t(k,j)==1) THEN
                       IF (qqm(i,j)>=qq(1)) THEN
                          DS(i,j)=0
                       ELSE IF (qqm(i,j)<qq(1)) THEN
                           DS(i,j)=1
                       END IF
                   END IF
                ELSE IF (levell(j)<-90 .and. levelc(j)>-90) THEN
                   kk=levelc(j)
                   IF (m_t(k,j)==-1) THEN
                       DS(i,j)=-99
                   ELSE IF (EXP(Ls(i,j,2))>=SUM(Cbar(1:kk)) .and. m_t(k,j)==1) THEN
                       
                       IF (qqm(i,j)>=qq(2)) THEN
                           DS(i,j)=1
                       ELSE IF (qqm(i,j)<qq(2)) THEN
                           DS(i,j)=0
                       END IF

                   ELSE IF (EXP(Ls(i,j,2))<SUM(Cbar(1:kk)) .and. m_t(k,j)==1) THEN
                       IF (qqm(i,j)>=qq(2)) THEN
                          DS(i,j)=0
                       ELSE IF (qqm(i,j)<qq(2)) THEN
                           DS(i,j)=1
                       END IF
                   END IF
                ELSE IF (levell(j)>-90 .and. levelc(j)>-90) THEN
                    kk=levelc(j)
                    tt=levell(j)
                   IF (m_t(k,j)==-1) THEN
                       DS(i,j)=-99
                   ELSE IF (EXP(Ls(i,j,1))>=SUM(Lbar(1:tt)) .and. EXP(Ls(i,j,2))>=SUM(Cbar(1:kk)) .and. m_t(k,j)==1) THEN
                       IF (qqm(i,j)>=MAX(qq(1),qq(2))) THEN
                           DS(i,j)=1
                       ELSE IF (qqm(i,j)<MAX(qq(1),qq(2))) THEN
                           DS(i,j)=0
                       END IF
                       
                   ELSE IF (EXP(Ls(i,j,1))<SUM(Lbar(1:tt)) .and. EXP(Ls(i,j,2))>=SUM(Cbar(1:kk)) .and. m_t(k,j)==1) THEN
                       
                       IF (qqm(i,j)>=MAX(qq(1),qq(2))) THEN
                           DS(i,j)=0
                       ELSE IF (qqm(i,j)<MAX(qq(1),qq(2))) THEN
                           DS(i,j)=1
                       END IF
                       
                   ELSE IF (EXP(Ls(i,j,1))>=SUM(Lbar(1:tt)) .and. EXP(Ls(i,j,2))<SUM(Cbar(1:kk)) .and. m_t(k,j)==1) THEN
                       
                       IF (qqm(i,j)>=MAX(qq(1),qq(2))) THEN
                           DS(i,j)=0
                       ELSE IF (qqm(i,j)<MAX(qq(1),qq(2))) THEN
                           DS(i,j)=1
                       END IF
                   END IF
                    
                ELSE IF (levell(j)<-90 .and. levelc(j)<-90 .and. levelf(j)>-90) THEN
                   kk=levelf(j)
                   IF (m_t(k,j)==-1) THEN
                       DS(i,j)=-99
                   ELSE IF (EXP(Ls(i,j,3))>=SUM(fbar(1:kk)) .and. m_t(k,j)==1) THEN
                       IF (qqm(i,j)>=qq(3)) THEN
                           DS(i,j)=1
                       ELSE IF (qqm(i,j)<qq(3)) THEN
                           DS(i,j)=0
                       END IF
                       
                   ELSE IF (EXP(Ls(i,j,3))<SUM(fbar(1:kk)) .and. m_t(k,j)==1) THEN
                       IF (qqm(i,j)>=qq(3)) THEN
                           DS(i,j)=0
                       ELSE IF (qqm(i,j)<qq(3)) THEN
                           DS(i,j)=1
                       END IF
                       
                   END IF
                END IF    
            END DO
        END DO
                
         !!!!!!!!!!!!!!!!!From here !!!!!!!!!!!!!!!!!!!!!!!
        !!!!!!!!!!!!!!!!!From here !!!!!!!!!!!!!!!!!!!!!!!
    !!! language tasks!!!!
    indext=0
    DO i=1,N_inds
      DO j=1,N_testt
         IF (ds(i,j)>-90 .and. levelg(j)<-90) THEN
             indext(j,1)=indext(j,1)+1
             IF (ds(i,j)>0.5) indext(j,2)=indext(j,2)+1
         END IF         
      END DO
    END DO
    
    DO i=1,N_testt
        IF (indext(i,1)>0.5) THEN
            mom_simu(i)=DBLE(indext(i,2))/DBLE(indext(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
        ELSE 
            mom_simu(i)=-99
        END IF
    END DO
    
    !!!! generate artl variables !!!!
    artls=-1
    DO i=1,N_inds
        DO j=1,N_testt
            IF (ds(i,j)>-90 .and. levell(j)>-90) THEN
                kk=levell(j)
                k=i-INT((i-1)/N_IND)*N_ind
                m=repl(k,j)
                artls(i,kk,m)=ds(i,j)
                inters_l(i,kk,m,1:3)=x_testt(k,j,1:3)
            END IF
        END DO
    END DO
    
    artcs=-1
    DO i=1,N_inds
        DO j=1,N_testt
            IF (ds(i,j)>-90 .and. levelc(j)>-90) THEN
                kk=levelc(j)
                k=i-INT((i-1)/N_IND)*N_ind
                m=repc(k,j)
                artcs(i,kk,m)=ds(i,j)
                inters_c(i,kk,m,1:3)=x_testt(k,j,1:3)
            END IF
        END DO
    END DO
    
    artfs=-1
    DO i=1,N_inds
        DO j=1,N_testt
            IF (ds(i,j)>-90 .and. levelf(j)>-90) THEN
                kk=levelf(j)
                k=i-INT((i-1)/N_IND)*N_ind
                m=repf(k,j)
                artfs(i,kk,m)=ds(i,j)
                inters_f(i,kk,m,1:3)=x_testt(k,j,1:3)
            END IF
        END DO
    END DO
    
    !artgs=-1
    !DO i=1,N_inds
    !    DO j=1,N_testt
    !        IF (ds(i,j)>-90 .and. levelg(j)>-90) THEN
    !            kk=levelg(j)
    !            k=i-INT((i-1)/N_IND)*N_ind
    !            m=repg(k,j)
    !            artgs(i,kk,m)=ds(i,j)
    !        END IF
    !    END DO
    !END DO
    
    !!!! passing rate at each level!!!!

    DO j=N_testt+1,N_testt+N_levell
       m=0
       n=0 
       DO i=1,N_inds
           DO k=1,20
               IF (artls(i,j-N_testt,k)>-0.5) THEN
                   m=m+1
                   IF (artls(i,j-N_testt,k)==1) n=n+1
               END IF
           END DO 
       END DO
       IF (m>0) THEN
           mom_simu(j)=DBLE(n)/DBLE(m)
       ELSE
           mom_simu(j)=-99
       END IF
       
    END DO
    
    !!!! correlation across levels !!!!
    
    DO j=2,N_levell-1
      DO t=1,5
        DO k=1,5
          m=0
          n=0 
          DO i=1,N_inds 
             IF (artls(i,j,t)==1 .and. artls(i,j+1,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artls(i,j,t)==1 .and. artls(i,j+1,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(N_testt+N_levell+(j-2)*25+(t-1)*5+k)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(N_testt+N_levell+(j-2)*25+(t-1)*5+k)=-99
          END IF
        END DO
      END DO
   END DO        
        
   DO j=2,N_levell
        l=0
      DO t=1,4
        DO k=t+1,5
          m=0
          n=0 
          l=l+1
          DO i=1,N_inds 
             IF (artls(i,j,t)==1 .and. artls(i,j,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artls(i,j,t)==1 .and. artls(i,j,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(25*(N_levell-2)+N_testt+N_levell+(j-2)*10+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(25*(N_levell-2)+N_testt+N_levell+(j-2)*10+l)=-99
          END IF
        END DO
      END DO
  END DO       
  
  DO j=2,N_levell-2
      DO t=1,5
        DO k=1,5
          m=0
          n=0 
          DO i=1,N_inds 
             IF (artls(i,j,t)==1 .and. artls(i,j+2,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artls(i,j,t)==1 .and. artls(i,j+2,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(25*(N_levell-2)+N_testt+N_levell+(N_levell-1)*10+(j-2)*25+(t-1)*5+k)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(25*(N_levell-2)+N_testt+N_levell+(N_levell-1)*10+(j-2)*25+(t-1)*5+k)=-99
          END IF
        END DO
      END DO
  END DO
  
  
      indext=0   ! task passing rate for new enrolled children
    DO i=1,N_inds
        k=i-INT((i-1)/N_IND)*N_ind
      DO j=1,N_testt
          as=DBLE(firstt(k,2)-month(j)-DBLE(time(j))/4.0d0)
         IF (ds(i,j)>-90 .and. levelg(j)<-90 .AND. as>-1.25) THEN
             indext(j,1)=indext(j,1)+1
             IF (ds(i,j)>0.5) indext(j,2)=indext(j,2)+1
         END IF         
      END DO
    END DO
    
    DO i=1,N_testt
        IF (indext(i,1)>0.5) THEN
            mom_simu(25*(2*N_levell-5)+N_testt+N_levell+(N_levell-1)*10+i)=DBLE(indext(i,2))/DBLE(indext(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
        ELSE 
            mom_simu(25*(2*N_levell-5)+N_testt+N_levell+(N_levell-1)*10+i)=-99
        END IF
    END DO
    
    
    indext=0   ! task passing rate for old enrolled children
    DO i=1,N_inds
        k=i-INT((i-1)/N_IND)*N_ind
      DO j=1,N_testt
          as=DBLE(firstt(k,2)-month(j)-DBLE(time(j))/4.0d0)
         IF (ds(i,j)>-90 .and. levelg(j)<-90 .AND. as<-1.25) THEN
             indext(j,1)=indext(j,1)+1
             IF (ds(i,j)>0.5) indext(j,2)=indext(j,2)+1
         END IF         
      END DO
    END DO
    
    DO i=1,N_testt
        IF (indext(i,1)>0.5) THEN
            mom_simu(25*(2*N_levell-5)+2*N_testt+N_levell+(N_levell-1)*10+i)=DBLE(indext(i,2))/DBLE(indext(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
        ELSE 
            mom_simu(25*(2*N_levell-5)+2*N_testt+N_levell+(N_levell-1)*10+i)=-99
        END IF
    END DO
    
    !!!! passing rate for cognitive tasks at each level!!!!
    
    kk=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10

    DO j=kk+1,kk+N_levelc
       m=0
       n=0 
       DO i=1,N_inds
           DO k=1,11
               IF (artcs(i,j-kk,k)>-0.5) THEN
                   m=m+1
                   IF (artcs(i,j-kk,k)==1) n=n+1
               END IF
           END DO 
       END DO
       IF (m>0) THEN
           mom_simu(j)=DBLE(n)/DBLE(m)
       ELSE
           mom_simu(j)=-99
       END IF
    END DO
    !!!! passing rate for fine motor tasks at each level!!!!
    
    kk=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc

    DO j=kk+1,kk+N_levelf
       m=0
       n=0 
       DO i=1,N_inds
           DO k=1,6
               IF (artfs(i,j-kk,k)>-0.5) THEN
                   m=m+1
                   IF (artfs(i,j-kk,k)==1) n=n+1
               END IF
           END DO 
       END DO
       IF (m>0) THEN
           mom_simu(j)=DBLE(n)/DBLE(m)
       ELSE
           mom_simu(j)=-99
       END IF
       
    END DO
    
    !    !!!! passing rate for gross motor tasks at each level!!!!
    !
    !kk=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf
    !
    !DO j=kk+1,kk+N_levelg
    !   m=0
    !   n=0 
    !   DO i=1,N_inds
    !       DO k=1,3
    !           IF (artgs(i,j-kk,k)>-0.5) THEN
    !               m=m+1
    !               IF (artgs(i,j-kk,k)==1) n=n+1
    !           END IF
    !       END DO 
    !   END DO
    !   IF (m>0) THEN
    !       mom_simu(j)=DBLE(n)/DBLE(m)
    !   ELSE
    !       mom_simu(j)=-99
    !   END IF
    !END DO
    
    
    
    !!!!! correlation across levels
    ! cognitive skills
    
    mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf
    
    l=0
    DO j=1,N_levelc-2
        IF (j<=2 .or. (j>4 .and. j<=8)) THEN
            tt=5
        ELSE IF (j>2 .and. j<5) THEN
            tt=2
        ELSE IF (j==9) THEN
            tt=2
        ELSE IF (j==10) THEN
            tt=4
        ELSE IF (j==11) THEN
            tt=5
        END IF
        
        IF (j<2 .or. (j>3 .and. j<8)) THEN
            kk=5
        ELSE IF (j>1 .and. j<4) THEN
            kk=2
        ELSE IF (j==8) THEN
            kk=2
        ELSE IF (j==9) THEN
            kk=4
        ELSE IF (j==10) THEN
            kk=5
        ELSE IF (j==11) THEN
            kk=3
        END IF
        
         DO t=1,tt
           DO k=1,kk
             m=0
             n=0 
             l=l+1
             DO i=1,N_inds 
                IF (artcs(i,j,t)==1 .and. artcs(i,j+1,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                    m=m+1
                    IF (artcs(i,j,t)==1 .and. artcs(i,j+1,k)==1) n=n+1
                END IF
             END DO
             
             IF (m>0) THEN
                mom_simu(mmm+l)=DBLE(n)/DBLE(m)
             ELSE
                mom_simu(mmm+l)=-99
             END IF
           END DO
         END DO
    END DO
    
    mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177
    
   l=0 
   DO j=1,N_levelf-1
       IF (j==1) THEN
           tt=5
           kk=2
       ELSE IF (j==2) THEN
           tt=2
           kk=5
       ELSE IF (j==3) THEN
           tt=5
           kk=5
       ELSE IF (j==4) THEN
           tt=5
           kk=3
       ELSE IF(j==5) THEN
           tt=3
           kk=4
       ELSE IF (j==6) THEN
           tt=4
           kk=3
       END IF
      DO t=1,tt
        DO k=1,kk
          m=0
          n=0
          l=l+1
          DO i=1,N_inds 
             IF (artfs(i,j,t)==1 .and. artfs(i,j+1,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artfs(i,j,t)==1 .and. artfs(i,j+1,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(mmm+l)=-99
          END IF
        END DO
      END DO
   END DO
   
   mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84
   
   !! gross motor skill
   !
   !l=0 
   !DO j=4,N_levelg-1
   !    IF (j==4) THEN
   !        tt=3
   !        kk=2
   !    ELSE IF (j==5) THEN
   !        tt=2
   !        kk=2
   !    ELSE IF (j==6) THEN
   !        tt=2
   !        kk=3
   !    ELSE IF (j==7) THEN
   !        tt=3
   !        kk=3
   !    ELSE IF(j==8) THEN
   !        tt=3
   !        kk=2
   !    END IF
   !   DO t=1,tt
   !     DO k=1,kk
   !       m=0
   !       n=0
   !       l=l+1
   !       DO i=1,N_inds 
   !          IF (artgs(i,j,t)==1 .and. artgs(i,j+1,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
   !              m=m+1
   !              IF (artgs(i,j,t)==1 .and. artgs(i,j+1,k)==1) n=n+1
   !          END IF
   !       END DO
   !       
   !       IF (m>0) THEN
   !          mom_simu(mm+l)=DBLE(n)/DBLE(m)
   !       ELSE
   !          mom_simu(mm+l)=-99
   !       END IF
   !     END DO
   !   END DO
   !END DO
   
   !mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84
   
    ! cognitive
    
    l=0
    DO j=1,N_levelc-1
        IF (j<=2 .or. (j>4 .and. j<=8)) THEN
            tt=4
            kk=5
        ELSE IF (j>2 .and. j<5) THEN
            tt=1
            kk=2
        ELSE IF (j==9) THEN
            tt=1
            kk=2
        ELSE IF (j==10) THEN
            tt=3
            kk=4
        ELSE IF (j==11) THEN
            tt=4
            kk=5
        ELSE IF (j==12) THEN
            tt=2
            kk=3
        END IF
        
         DO t=1,tt
           DO k=t+1,kk
             m=0
             n=0 
             l=l+1
             DO i=1,N_inds 
                IF (artcs(i,j,t)==1 .and. artcs(i,j,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                    m=m+1
                    IF (artcs(i,j,t)==1 .and. artcs(i,j,k)==1) n=n+1
                END IF
             END DO
             
             IF (m>0) THEN
                mom_simu(mmm+l)=DBLE(n)/DBLE(m)
             ELSE
                mom_simu(mmm+l)=-99
             END IF
           END DO
         END DO
    END DO
    
    mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82
    
    ! fine motor
    
    l=0
    DO j=1,N_levelf
        
       IF (j==1) THEN
           tt=4
           kk=5
       ELSE IF (j==2) THEN
           tt=1
           kk=2
       ELSE IF (j==3) THEN
           tt=4
           kk=5
       ELSE IF (j==4) THEN
           tt=4
           kk=5
       ELSE IF(j==5) THEN
           tt=2
           kk=3
       ELSE IF (j==6) THEN
           tt=3
           kk=4
       ELSE IF (j==7) THEN
           tt=2
           kk=3
       END IF
        
      DO t=1,tt
        DO k=t+1,kk
          m=0
          n=0
          l=l+1
          DO i=1,N_inds 
             IF (artfs(i,j,t)==1 .and. artfs(i,j,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artfs(i,j,t)==1 .and. artfs(i,j,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(mmm+l)=-99
          END IF
        END DO
      END DO
    END DO
    
    mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43
    
    !! gross motor
    !
    !l=0
    !DO j=4,N_levelg
    !    
    !   IF (j==4) THEN
    !       tt=2
    !       kk=3
    !   ELSE IF (j==5) THEN
    !       tt=1
    !       kk=2
    !   ELSE IF (j==6) THEN
    !       tt=1
    !       kk=2
    !   ELSE IF (j==7) THEN
    !       tt=2
    !       kk=3
    !   ELSE IF(j==8) THEN
    !       tt=2
    !       kk=3
    !   ELSE IF (j==9) THEN
    !       tt=1
    !       kk=2
    !   END IF
    !    
    !  DO t=1,tt
    !    DO k=t+1,kk
    !      m=0
    !      n=0
    !      l=l+1
    !      DO i=1,N_inds 
    !         IF (artgs(i,j,t)==1 .and. artgs(i,j,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
    !             m=m+1
    !             IF (artgs(i,j,t)==1 .and. artgs(i,j,k)==1) n=n+1
    !         END IF
    !      END DO
    !      
    !      IF (m>0) THEN
    !         mom_simu(mm+l)=DBLE(n)/DBLE(m)
    !      ELSE
    !         mom_simu(mm+l)=-99
    !      END IF
    !    END DO
    !  END DO
    !END DO
    
    !mm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+84+43
    
    ! cognitive
    
    l=0
    DO j=1,N_levelc-3
        IF (j==1) THEN
            tt=5
            kk=2
        ELSE IF (j==2) THEN
            tt=5
            kk=2
        ELSE IF (j==3) THEN
            tt=2
            kk=5
        ELSE IF (j==4) THEN
            tt=2
            kk=5
        ELSE IF (j==5) THEN
            tt=5
            kk=5
        ELSE IF (j==6) THEN
            tt=5
            kk=5
        ELSE IF (j==7) THEN
            tt=5
            kk=2
        ELSE IF (j==8) THEN
            tt=5
            kk=4
        ELSE IF (j==9) THEN
            tt=2
            kk=5
        ELSE IF (j==10) THEN
            tt=4
            kk=3
        END IF
        
        
      DO t=1,tt
        DO k=1,kk
          m=0
          n=0 
          l=l+1
          DO i=1,N_inds 
             IF (artcs(i,j,t)==1 .and. artcs(i,j+2,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artcs(i,j,t)==1 .and. artcs(i,j+2,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(mmm+l)=-99
          END IF
        END DO
      END DO
    END DO  
    
    mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43+142
    
     ! fine motor
    
    l=0
    DO j=1,N_levelf-2
        IF (j==1) THEN
            tt=5
            kk=5
        ELSE IF (j==2) THEN
            tt=2
            kk=5
        ELSE IF (j==3) THEN
            tt=5
            kk=3
        ELSE IF (j==4) THEN
            tt=5
            kk=4
        ELSE IF (j==5) THEN
            tt=3
            kk=3
        END IF
        
      DO t=1,tt
        DO k=1,kk
          m=0
          n=0 
          l=l+1
          DO i=1,N_inds 
             IF (artfs(i,j,t)==1 .and. artfs(i,j+2,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artfs(i,j,t)==1 .and. artfs(i,j+2,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(mmm+l)=-99
          END IF
        END DO
      END DO
    END DO  
    
    mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43+142+79
    
    !! gross motor
    !
    !l=0
    !DO j=4,N_levelg-2
    !    IF (j==4) THEN
    !        tt=3
    !        kk=2
    !    ELSE IF (j==5) THEN
    !        tt=2
    !        kk=3
    !    ELSE IF (j==6) THEN
    !        tt=2
    !        kk=3
    !    ELSE IF (j==7) THEN
    !        tt=3
    !        kk=2
    !    END IF
    !    
    !  DO t=1,tt
    !    DO k=1,kk
    !      m=0
    !      n=0 
    !      l=l+1
    !      DO i=1,N_inds 
    !         IF (artgs(i,j,t)==1 .and. artgs(i,j+2,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
    !             m=m+1
    !             IF (artgs(i,j,t)==1 .and. artgs(i,j+2,k)==1) n=n+1
    !         END IF
    !      END DO
    !      
    !      IF (m>0) THEN
    !         mom_simu(mm+l)=DBLE(n)/DBLE(m)
    !      ELSE
    !         mom_simu(mm+l)=-99
    !      END IF
    !    END DO
    !  END DO
    !END DO  
    
     !mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+84+43+142+79
 
    DULs_l=-99
    inds_l=-99
    DO i=1,N_inds
        DO k=1,N_levell
            DO j=1,20
                IF (artls(i,k,j)==1 .and. j==1) THEN
                    duls_l(i,k)=1
                    inds_l(i,k)=1
                ELSE IF (artls(i,k,j)==0 .and. j==1) THEN
                    duls_l(i,k)=j
                    inds_l(i,k)=0
                ELSE IF (j>1 ) THEN
                    IF( artls(i,k,j-1)==0 .and. artls(i,k,j)==0 .and. inds_l(i,k)==0) THEN
                       duls_l(i,k)=j
                       inds_l(i,k)=0
                       IF (j<20 .and. artls(i,k,j+1)==-1) THEN
                           duls_l(i,k)=j+1
                           inds_l(i,k)=1
                       END IF  
                    ELSE IF (j>1 .and. artls(i,k,j-1)==0 .and. artls(i,k,j)==0 .and. inds_l(i,k)==1) THEN
                        duls_l(i,k)=duls_l(i,k)
                        inds_l(i,k)=1
                    ELSE IF (j>1 .and. artls(i,k,j-1)==0 .and. artls(i,k,j)==1 .and. inds_l(i,k)==0) THEN
                        duls_l(i,k)=j
                        inds_l(i,k)=1
                    ELSE IF (j>1 .and. artls(i,k,j-1)==0 .and. artls(i,k,j)==1 .and. inds_l(i,k)==1) THEN
                        duls_l(i,k)=duls_l(i,k)
                        inds_l(i,k)=1   
                    END IF
                END IF
                
            END DO           
                
        END DO
    END DO
    
    DO k=1,n_levell
        m=0
        n=0
        DO i=1,n_inds
            IF (duls_l(i,k)>-90) THEN
               m=m+1
               n=n+duls_l(i,k)
            END IF
        END DO
        IF (m>0) THEN
          mom_simu(mmm+k)=DBLE(n)/DBLE(m)
        ELSE 
           mom_simu(mmm+k)=-99
        END IF
    END DO
    
    mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43+142+79+n_levell
    
    dus_l=0
    DO k=1,n_levell
        DO i=1,N_inds
            IF (duls_l(i,k)>-90) THEN
                dus_l(k)=dus_l(k)+1
            END IF
        END DO
    END DO
    DO k=1,n_levell
        allocate(avse(dus_l(k)))
        m=0
        DO i=1,n_inds
            IF (duls_l(i,k)>-90) THEN
                m=m+1
                avse(m)=duls_l(i,k)
            END IF
        END DO
        mom_simu(mmm+k)=STDEV_V(AVSE)
        DEALLOCATE(avse)
    END DO
     mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43+142+79+2*n_levell

    INDs_l=0
    interls_l=0.0d0
    DO k=2,N_levell
        DO i=1,N_inds
          DO j=1,20
              IF (inters_l(i,k,j,1)>-90) THEN
                  inds_l(i,k)=inds_l(i,k)+1
                  interls_l(i,k,1)=interls_l(i,k,1)+inters_l(i,k,j,1)
                  interls_l(i,k,2)=interls_l(i,k,2)+inters_l(i,k,j,2)
                  interls_l(i,k,3)=interls_l(i,k,3)+inters_l(i,k,j,3)
              END IF
          END DO
          IF (inds_l(i,k)>0) THEN
               interls_l(i,k,1)=DBLE(interls_l(i,k,1))/DBLE(inds_l(i,k))
               interls_l(i,k,2)=DBLE(interls_l(i,k,2))/DBLE(inds_l(i,k))
               interls_l(i,k,3)=DBLE(interls_l(i,k,3))/DBLE(inds_l(i,k))
          END IF
        END DO
    END DO
    
    dus_l=0
    DO k=1,n_levell
        DO i=1,N_inds
            IF (inds_l(i,k)>0 .AND. duls_l(i,k)>-90) THEN
            dus_l(k)=dus_l(k)+1          
            END IF
        END DO
    END DO
    
    DO k=1,n_levell
        allocate(cors(dus_l(k),4))
        m=0
        DO i=1,n_inds
             IF (inds_l(i,k)>0 .AND. duls_l(i,k)>-90) THEN
                m=m+1
                cors(m,1)=interls_l(i,k,1)
                cors(m,2)=interls_l(i,k,2)
                cors(m,3)=interls_l(i,k,3)
                cors(m,4)=duls_l(i,k)
             END IF
        END DO
        !if (k==2) write(*,*) m,Standard_Deviation(cors(:,4))
        IF (k>1 .and. Standard_Deviation(cors(:,4))>1.0d-6) THEN
        mom_simu(mmm+3*(k-1)+1)=Correlation(cors(:,1),cors(:,4))
        mom_simu(mmm+3*(k-1)+2)=Correlation(cors(:,2),cors(:,4))
        mom_simu(mmm+3*(k-1)+3)=Correlation(cors(:,3),cors(:,4))
        ELSE IF (k>1 .and. Standard_Deviation(cors(:,4))<1.0d-6) THEN
            mom_simu(mmm+3*(k-1)+1)=-3
            mom_simu(mmm+3*(k-1)+2)=-3
            mom_simu(mmm+3*(k-1)+3)=-3
        ELSE
            mom_simu(mmm+3*(k-1)+1)=-99
            mom_simu(mmm+3*(k-1)+2)=-99
            mom_simu(mmm+3*(k-1)+3)=-99
        END IF
        
            
        DEALLOCATE(cors)
    END DO
    mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43+142+79+5*n_levell
    
    
    !!! cognitive 
    
    
    DUcs_c=-99
    inds_c=-99
    DO i=1,N_inds
        DO k=1,N_levelc
            DO j=1,11
                IF (artcs(i,k,j)==1 .and. j==1) THEN
                    ducs_c(i,k)=1
                    inds_c(i,k)=1
                ELSE IF (artcs(i,k,j)==0 .and. j==1) THEN
                    ducs_c(i,k)=j
                    inds_c(i,k)=0
                ELSE IF (j>1 ) THEN
                    IF( artcs(i,k,j-1)==0 .and. artcs(i,k,j)==0 .and. inds_c(i,k)==0) THEN
                       ducs_c(i,k)=j
                       inds_c(i,k)=0
                       IF (j<11 ) THEN
                           IF (artcs(i,k,j+1)==-1) THEN
                           ducs_c(i,k)=j+1
                           inds_c(i,k)=1
                           END IF
                       END IF  
                    ELSE IF (j>1 .and. artcs(i,k,j-1)==0 .and. artcs(i,k,j)==0 .and. inds_c(i,k)==1) THEN
                        ducs_c(i,k)=ducs_c(i,k)
                        inds_c(i,k)=1
                    ELSE IF (j>1 .and. artcs(i,k,j-1)==0 .and. artcs(i,k,j)==1 .and. inds_c(i,k)==0) THEN
                        ducs_c(i,k)=j
                        inds_c(i,k)=1
                    ELSE IF (j>1 .and. artcs(i,k,j-1)==0 .and. artcs(i,k,j)==1 .and. inds_c(i,k)==1) THEN
                        ducs_c(i,k)=ducs_c(i,k)
                        inds_c(i,k)=1   
                    END IF
                END IF
                
            END DO           
                
        END DO
    END DO
    
    DO k=1,n_levelc
        m=0
        n=0
        DO i=1,n_inds
            IF (ducs_c(i,k)>-90) THEN
               m=m+1
               n=n+ducs_c(i,k)
            END IF
        END DO
        IF (m>0) THEN
          mom_simu(mmm+k)=DBLE(n)/DBLE(m)
        ELSE 
           mom_simu(mmm+k)=-99
        END IF
    END DO
    
    mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43+142+79+5*n_levell+n_levelc
    
    dus_c=0
    DO k=1,n_levelc
        DO i=1,N_inds
            IF (ducs_c(i,k)>-90) THEN
                dus_c(k)=dus_c(k)+1
            END IF
        END DO
    END DO
    DO k=1,n_levelc
        allocate(avse(dus_c(k)))
        m=0
        DO i=1,n_inds
            IF (ducs_c(i,k)>-90) THEN
                m=m+1
                avse(m)=ducs_c(i,k)
            END IF
        END DO
        mom_simu(mmm+k)=STDEV_V(AVSE)
        DEALLOCATE(avse)
    END DO
     mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43+142+79+5*n_levell+2*n_levelc

    INDs_c=0
    intercs_c=0.0d0
    DO k=1,N_levelc
        DO i=1,N_inds
          DO j=1,11
              IF (inters_c(i,k,j,1)>-90) THEN
                  inds_c(i,k)=inds_c(i,k)+1
                  intercs_c(i,k,1)=intercs_c(i,k,1)+inters_c(i,k,j,1)
                  intercs_c(i,k,2)=intercs_c(i,k,2)+inters_c(i,k,j,2)
                  intercs_c(i,k,3)=intercs_c(i,k,3)+inters_c(i,k,j,3)
              END IF
          END DO
          IF (inds_c(i,k)>0) THEN
               intercs_c(i,k,1)=DBLE(intercs_c(i,k,1))/DBLE(inds_c(i,k))
               intercs_c(i,k,2)=DBLE(intercs_c(i,k,2))/DBLE(inds_c(i,k))
               intercs_c(i,k,3)=DBLE(intercs_c(i,k,3))/DBLE(inds_c(i,k))
          END IF
        END DO
    END DO
    
    dus_c=0
    DO k=1,n_levelc
        DO i=1,N_inds
            IF (inds_c(i,k)>0 .AND. ducs_c(i,k)>-90) THEN
            dus_c(k)=dus_c(k)+1          
            END IF
        END DO
    END DO
    
    DO k=1,n_levelc
        allocate(cors(dus_c(k),4))
        m=0
        DO i=1,n_inds
             IF (inds_c(i,k)>0 .AND. ducs_c(i,k)>-90) THEN
                m=m+1
                cors(m,1)=intercs_c(i,k,1)
                cors(m,2)=intercs_c(i,k,2)
                cors(m,3)=intercs_c(i,k,3)
                cors(m,4)=ducs_c(i,k)
             END IF
        END DO
        IF (k<n_levelc) THEN
        mom_simu(mmm+3*(k-1)+1)=Correlation(cors(:,1),cors(:,4))
        mom_simu(mmm+3*(k-1)+2)=Correlation(cors(:,2),cors(:,4))
        mom_simu(mmm+3*(k-1)+3)=Correlation(cors(:,3),cors(:,4))
        ELSE 
            mom_simu(mmm+3*(k-1)+1)=-99
            mom_simu(mmm+3*(k-1)+2)=-99
            mom_simu(mmm+3*(k-1)+3)=-99
        END IF
        
            
        DEALLOCATE(cors)
    END DO

    mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43+142+79+5*n_levell+5*n_levelc   
     
     
     
    DUfs_f=-99
    inds_f=-99
    DO i=1,N_inds
        DO k=1,N_levelf
            DO j=1,6
                IF (artfs(i,k,j)==1 .and. j==1) THEN
                    dufs_f(i,k)=1
                    inds_f(i,k)=1
                ELSE IF (artfs(i,k,j)==0 .and. j==1) THEN
                    dufs_f(i,k)=j
                    inds_f(i,k)=0
                ELSE IF (j>1 ) THEN
                    IF( artfs(i,k,j-1)==0 .and. artfs(i,k,j)==0 .and. inds_f(i,k)==0) THEN
                       dufs_f(i,k)=j
                       inds_f(i,k)=0
                       IF (j<6 ) THEN
                           IF (artfs(i,k,j+1)==-1) THEN
                           dufs_f(i,k)=j+1
                           inds_f(i,k)=1
                           END IF
                       END IF  
                    ELSE IF (j>1 .and. artfs(i,k,j-1)==0 .and. artfs(i,k,j)==0 .and. inds_f(i,k)==1) THEN
                        dufs_f(i,k)=dufs_f(i,k)
                        inds_f(i,k)=1
                    ELSE IF (j>1 .and. artfs(i,k,j-1)==0 .and. artfs(i,k,j)==1 .and. inds_f(i,k)==0) THEN
                        dufs_f(i,k)=j
                        inds_f(i,k)=1
                    ELSE IF (j>1 .and. artfs(i,k,j-1)==0 .and. artfs(i,k,j)==1 .and. inds_f(i,k)==1) THEN
                        dufs_f(i,k)=dufs_f(i,k)
                        inds_f(i,k)=1   
                    END IF
                END IF
                
            END DO           
                
        END DO
    END DO
    
    DO k=1,n_levelf
        m=0
        n=0
        DO i=1,n_inds
            IF (dufs_f(i,k)>-90) THEN
               m=m+1
               n=n+dufs_f(i,k)
            END IF
        END DO
        IF (m>0) THEN
          mom_simu(mmm+k)=DBLE(n)/DBLE(m)
        ELSE 
           mom_simu(mmm+k)=-99
        END IF
    END DO
    
    mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43+142+79+5*n_levell+5*n_levelc +n_levelf  
    
        dus_f=0
    DO k=1,n_levelf
        DO i=1,N_inds
            IF (dufs_f(i,k)>-90) THEN
                dus_f(k)=dus_f(k)+1
            END IF
        END DO
    END DO
    DO k=1,n_levelf
        allocate(avse(dus_f(k)))
        m=0
        DO i=1,n_inds
            IF (dufs_f(i,k)>-90) THEN
                m=m+1
                avse(m)=dufs_f(i,k)
            END IF
        END DO
        mom_simu(mmm+k)=STDEV_V(AVSE)
        DEALLOCATE(avse)
    END DO
     mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43+142+79+5*n_levell+5*n_levelc +2*n_levelf  

    INDs_f=0
    interfs_f=0.0d0
    DO k=1,N_levelf
        DO i=1,N_inds
          DO j=1,6
              IF (inters_f(i,k,j,1)>-90) THEN
                  inds_f(i,k)=inds_f(i,k)+1
                  interfs_f(i,k,1)=interfs_f(i,k,1)+inters_f(i,k,j,1)
                  interfs_f(i,k,2)=interfs_f(i,k,2)+inters_f(i,k,j,2)
                  interfs_f(i,k,3)=interfs_f(i,k,3)+inters_f(i,k,j,3)
              END IF
          END DO
          IF (inds_f(i,k)>0) THEN
               interfs_f(i,k,1)=DBLE(interfs_f(i,k,1))/DBLE(inds_f(i,k))
               interfs_f(i,k,2)=DBLE(interfs_f(i,k,2))/DBLE(inds_f(i,k))
               interfs_f(i,k,3)=DBLE(interfs_f(i,k,3))/DBLE(inds_f(i,k))
          END IF
        END DO
    END DO
    
    dus_f=0
    DO k=1,n_levelf
        DO i=1,N_inds
            IF (inds_f(i,k)>0 .AND. dufs_f(i,k)>-90) THEN
            dus_f(k)=dus_f(k)+1          
            END IF
        END DO
    END DO
    
    DO k=1,n_levelf
        allocate(cors(dus_f(k),4))
        m=0
        DO i=1,n_inds
             IF (inds_f(i,k)>0 .AND. dufs_f(i,k)>-90) THEN
                m=m+1
                cors(m,1)=interfs_f(i,k,1)
                cors(m,2)=interfs_f(i,k,2)
                cors(m,3)=interfs_f(i,k,3)
                cors(m,4)=dufs_f(i,k)
             END IF
        END DO
        !IF (k>1) THEN
        mom_simu(mmm+3*(k-1)+1)=Correlation(cors(:,1),cors(:,4))
        mom_simu(mmm+3*(k-1)+2)=Correlation(cors(:,2),cors(:,4))
        mom_simu(mmm+3*(k-1)+3)=Correlation(cors(:,3),cors(:,4))
        !ELSE 
        !    mom_simu(mmm+3*(k-1)+1)=-99
        !    mom_simu(mmm+3*(k-1)+2)=-99
        !    mom_simu(mmm+3*(k-1)+3)=-99
        !END IF
        
            
        DEALLOCATE(cors)
    END DO
 
        mmm=25*(2*N_levell-5)+3*N_testt+N_levell+(N_levell-1)*10+n_levelc+n_levelf+177+84+82+43+142+79+5*n_levell+5*n_levelc +5*n_levelf 
    
    
    
    
    
    DO i=1,N_mom
        IF (ABS(moment(i))<1.0d-10) THEN
            dif(i,1)=((mom_simu(i)-moment(i)))**2.0D0
        ELSE 
            dif(i,1)=((DBLE((mom_simu(i)-moment(i))))**2.0D0)
        END IF 
   
    END DO
    
    !mm=MATMUL(MATMUL(TRANSPOSE(DIF),WEI),DIF)
    mom=SUM(dif)
    
    
    
                OPEN(6424,file="dif.out")
                  DO i = 1,N_mom
                    WRITE(6424,'(200F16.8)') SQRT(dif(i,1)),mom_simu(i),moment(i)
                  END DO
                CLOSE(6424)
                
   
   END FUNCTION mom
    
       
    
END MODULE smm